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Abstract 


In papers I and II we have presented the results of the most updated 1 -point closure 

model for the turbulent vertical diffusivities of momentum, heat and salt, K m ^ g . In this 

paper, In this paper, we derive the analytic expressions for K m ^ g using a new 2-point 

closure model that has recently been developed and successfully tested against some ~80 

turbulence statistics for different flows. The new model has no free parameters. The 

expressions for K , are analytical functions of two stability parameters: the Turner 
m,n,b 

number R^ (salinity gradient/temperature gradient) and the Richardson number Ri 
(temperature gradient /shear). The turbulent kinetic energy K and its rate of dissipation 
may be taken local or non-local (K-e model). 

Contrary to all previous models that to describe turbulent mixing below the mixed 
layer (ML) have adopted three adjustable "background diffusivities" for momentum, heat 
and salt, we propose a model that avoids such adjustable diffusivities. We assume that 
below the ML, K m ^ g have the same Junctional dependence on Ri and R ^ derived from the 
turbulence model. However, in order to compute Ri below the ML, we use data of vertical 
shear due to wave— breaking measured by Gargett et al. (1981). The procedure frees the 
model from adjustable background diffusivities and indeed we use the same model 
throughout the entire vertical extent of the ocean. 

Using the new K , , we run an O-GCM and present a variety of results that we 

m,n,s 

compare with Levitus and the KPP model. 

Since the traditional 1— point (used in papers I and II) and the new 2— point closure 
models used here represent different modeling philosophies and procedures, testing them in 
an O-GCM is indispensable. The basic motivation is to show that the new 2-point closure 
model gives results that are overall superior to the 1— point closure in spite of the fact that 
the latter rely on several adjustable parameters while the new 2— point closure has none. 

After the extensive comparisons presented in papers I and II, we conclude that the 
new model presented here is overall superior for it not only is parameter free but also 
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because is part of a more general turbulence model that has been previously successfully 
tested on a wide variety of other types of turbulent flows. 
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I.Introduction 


a) O-GCM dynamic equations. Turbulent Diffusivities 

The dynamic equations for the large scale mean velocity, temperature and salinity, 
Ik, T and S, are given by (Semtner; 1995; Semtner and Chervin, 1992; Anderson and 
Willebrand, 1992) 


D tt , d _ dP . 
Dt u i + ^. r ij- _ 3x + -• 

DT , d . _ d , JTY , ^ 
Dt + dx“\ ~ + - 

DS d _ d ( as X 
Dt + c?x. s i - + - 


(la) 

(lb) 

(lc) 


where D/Dt=d/dt+\J^d/dXy \ and « are the temperature and salt kinematic diffusivities. 


Eqs.(la— c) depend on turbulence via the the Reynolds stresses r^UjUj and the heat and 
salt fluxes tn=u79, Sj=Uj<7, where u-, 6 and a represent the turbulent (fluctuating) 
components of the velocity, temperature and salinity fields. The dots represent external 
forces, rotation, etc. The vertical turbulent diffusivities K , (cm 2 s 4 ) enter via the 

relations: 


— i/ 3U —a „ 91 — u dS n .v 

uw = - K m?F' Wd= ~ K hW V,c = - K sm (ld) 

Historically, a value of K^=l was suggested by Munk (1966). Bryan (1991) and Colin de 

2 j 

Verdiere (1988) later showed that the polar heat transport scales like K^ /3 and that in 
order to reproduce its observed value (~lPetawatt), 0— GCM require K^~2 (Holland, 1989; 
Mellor, 1989; Cummings et., 1990). Chen et al. (1994) concluded that an improved mixed 
layer ameliorates the large scale equatorial circulation as well as the annual cycles of the 
SST (sea surface temperatures) while on a global scale McWilliams (1996) has recently 
stressed that O-GCMs are sensitive, among other things, to the vertical diffusivities. From 
high resolution profiles and tracer experiments (Mourn 1989; Gargett, 1989; Ledwell et al. 
1993; McPhee and Martison 1994; Toole et al. 1994; Polzin et al. 1997), it was concluded 
that in the ocean mixed layer can be as large as ( 1 — 2) 10 2 , while in the thermocline 


4 



Kj i ~(0. 11-0.6) with perhaps the tendency to vary as N" 1 , where N is the Brunt-Vaisala 
frequency (Gargett and Holloway 1984; Gargett et ah, 1989; Gargett 1984, 1993; Matear 
and Wong, 1997). Thus, diffusivities are not constant with depth and different processes 
contribute to the K's at different depths resulting in different values. 


b) general constraints on the K’s 

Since the dimensions of the diffusivities are (length) 2 (time)' 1 , to construct the 
diffusivities K's one needs two independent variables. Over the years, three expressions 
have been adopted: 


1) Richardson law: 

Diffusivity ~ t ^ 3 f ^ 3 

(le) 

2) Prandtl law: 

Diffusivity ~ YSt 

(If) 

3) Turbulence modeling: 

r2 

Diffusivity ~ — 

(lg) 


where i is a typical length scale, K is the turbulent kinetic energy and e is the rate of 
dissipation of K. Richardson law (a precursor of Kolmogorov law) emphasizes the diffusive 
nature of large eddies while Prandtl law is reminiscent of kinetic theory of gases and mean 
free path arguments. Both expressions entail a length scale l which is usually difficult to 
quantify, especially in the case of stable stratification (Cheng and Canuto, 1994), whereas 
(lg) depends on K and e for which we shall derive two differential equations thus avoiding 
the need of an l. On that basis alone, (lg) is preferable to both (le— f). Even so, (lg) is still 
not satisfactory. In fact, consider the case of stable stratification when shear is the source of 


turbulence. In the absence of salt, the turbulent Prandtl number, 

»T = Kf < lh > 

is known to increase with stratification represented by the Richardson number Ri, 

a T = 3Ri <7 T (R'i) > 0 (li) 

(Webster, 1964; Istweire and Helland, 1989; Schumann, 1991; Wang et al., 1996). If we 


limit the functional dependence of the K's to the forms (le-g), we obtain <r T =constant 
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which violates the second of (li). For example, the model of Pacanoswki and Philander 
does not satisfy (li) since it yields a constant a . The recent model for K , of D'Alessio 

1 1 II ^ XI 

et al. (1998) does not satisfy (li) either. In addition, the function S m must satisfy the well 
known asymptotic limits of RhO (shear driven flows) 

Ri->0: 2S m «0.1 (lj) 

as studies of shear flow indicate (Launder et al., 1975). 

Finally, the model must provide a correct value of the critical Richardson numbe: 
Ri cr above which turbulence can no longer be sustained by shear. While linear stability 
analysis implies that turbulence exists only for 

Ri< Ri a , Ri cr = 1/4 (lk) 

even early laboratory data by Taylor (cited in Monin and Yaglom, 1971, Vol.I, pages 
502—503) show that considerable turbulent exchange exists up to Ri~10. More recently, 
Martin (1985), Smart (1988) and Kundu et al. (1991) have found considerable mixing in 
the ocean up to Ri~l, a result also validated by the LES results of Wang et al. (1996). This 
means that a model for stably stratified turbulence must allow turbulence past Ri cr =l/4. 
The original MY model (Mellor and Yamada, 1982) does not satisfy this criterion since it 
yields Ri cr =0.19 nor do any of the subsequent models that have not significantly improved 
the physical content of the MY model. This brief discussion about a T (Ri) and Ri cr suffices 
to show that there are some basic data that any model must satisfy before being considered 
a candidate for ocean turbulence. 


The general functional form of the turbulent diffusivities for momentum, heat and salt 


IC u „ are of the form: 
m,h,s 


K = 2 li 2 S 
*m,h,s L e m,h,s 


where K and t are the turbulent kinetic energy and its rate of dissipation which in principle 
are given by two dynamic equations (the K— f model). The dimensionless structure 
Junctions S m ^ g must differ from one another so that: 
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In general we can write: 


W K S 


S m,h, s = W™’ «r"‘. “s VS > 


where a and a are the volume expansion coefficients a —-p A dp/ cfT and a =p A dp/ d S and 
where the shear VU can be generated either by external sources like in the mixed layer ML 
or by internal wave-breaking processes below the ML. If one introduces the Turner number 
and the Richardson number Ri: 


R p ~ S a s |f Ri “ N h /N u 


where 


N h = 


N 2 = 2(E..E..) 
u v ij ij ; 


N2 -~p^~ ~ S a sll " N h (1_R /> ) 


we can rewrite (2c) more concisely as 

S , = S , (R , Ri) (2f) 

m,h,s m,h,s v p ' v ' 

(Clearly, we could have also defined Ri in terms of N 2 rather than just the thermal 
gradient. We have chosen the latter for reasons of presentation of the results). One must 
distinguish the following four cases: 

SF ( salt fingers, salty -warm over fresh- cold): 

dS >n 3T >n 

3F >0 ’ 3F >0 > 

R >0, Ri>0 
P 

R^<1 Stable, N 2 >0, R >1 Unstable, N 2 <0 (3a) 

DC ( diffusive convection, fresh-cold over salty -warm): 

dS <(] dT <(] 

m <0 ’ W <0 ' 

R. >0, Ri<0 
P 

R <1 Unstable, N 2 <0, R >1 Stable, N 2 >0 (3b) 

r r 

DS ( doubly stable, fresh-warm over salty-cold) 

d$ <() dr 0 
3i <0 ’ ^->0, 

R^<0, Ri>0, N 2 >0, Stable (3c) 
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DU ( doubly unstable . salty-cold over fresh- warm): 


3S n 0T <(} 

^> 0 , -^<U, 

R^<0, Ri<0, N 2 <0, Unstable 


(3d) 


The stability /instability is predicated on the Brunt— Vaisala frequency N with N 2 >0 
(stable) and N 2 <0 (unstable). 

The general problem is to construct (2f) so as to encompass all four cases /3a -d/ 
First, there is ample evidence from laboratory and oceanic field data that show that K^ is 
different from K . In the SF case, the ratio K /K, >1 (Hamilton et al., 1989) while in the 
DC case, K h /K g >l (Kelley, 1984). Schmitt (1981) has also shown that the observed T-S 
relationship is not consistent with K^=K g . For a discussion and review of the importance 
of these processes and their extent in different parts of the ocean, see Turner (1967, 1973, 
1985), Schmitt (1994) and Zhang et al. (1998). In spite of this evidence, almost all O-GCM 
still assume 


K s =K h (3e) 

Recently, attempts have been made to overcome (3e) but the task is not easy. The main 
difficulty is that in the absence of a model capable of encompassing all four cases, SF, DC. 
DS and DU, the only alternative is to employ laboratory and ocean data to build the 
functional form of the diffusivities to be used in an 0— GCM. This is the approach 
employed by Large et al. (1994), Zhang et al. (1998, ZSH) and Merryfield et al. (1999, 
MHG) who used relations by Schmitt (1981) and Kelley (1984, 1990), among others. 

There is, however, an internal limitation to such a procedure since the available data 
refer to SF and DC but not to DS and DU which are also important (Duffy and Caldera, 
1999). Thus, away from the regions where SF and DC are active, the above authors take 

K m,h, S < DS ’ DlJ ) =° (3f) 

or, more precisely, they use a background diffusivity which is chosen primarily on grounds 

of numerical stability but whose physical origin must be an internal-wave breaking 

phenomenon. This is clearly not a satisfactory situation especially in view of the fact that 


8 



Since the studies by ZHS and MHG have shown the importance of double diffusion, the 
above procedure is certainly better than (3e) but still not fully satisfactory. The goal of this 
paper is to consider the same problem but with a different methodology. 

We develop a turbulence model to compute the three diffusivities for momentum, heat 
and salt to encompass the four processes SF, DC, DS and DU in the presence of an 
arbitrary shear. The inclusion of shear is quite relevant since is known to hamper the SF 
mechanism (Linden, 1971, 1974a, b; Kunze, 1990) and yet, the above procedures do not 
account for shear since they expressed and K g in terms of only one stability parameter 

R , rather than R and the Richardson number Ri. 

P P 

The model comes in three forms: 1) K and e are solutions of two dynamical equations, 
the so called K — t model whose main but not unique advantage is to avoid the introduction 
of any mixing length, 2) only one of them satisfies a differential equation while the other is 
taken to be the local limit of its dynamic equation and 3) both K and f are taken as the 
local limit of their respective dynamic equations. As we shall show, in model 3) all the 
relations are algebraic and one must solve a cubic equation. All the numerical results 
correspond to case 3). 

c) How to compute the diffusivities K’s 

Since from the first 0— GCM it clearly emerged that several key ocean features, e.g., 
polar heat transport, meridian circulation, T-S relationship (temperature— salinity), etc., 
depend sensitively on the vertical turbulent diffusivities of momentum K , heat and 
salt K (F.O. Bryan, 1987; K. Bryan, 1991) even though until recently all codes employed 
the approximation K g =K^. The need to reliably compute such diffusivities has been dealt 
in the past with three quite distinct approaches. 

d) descriptive-diagnostic method 

One can adopt a descriptive— diagnostic approach whereby the 0— GCM are used to 
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show that in order to obtain overall agreement with the data, the diffusivities K's must lie 
in a restricted range of values. Such a diagnostic, a posteriori approach is however 
inadequate in climate studies where the predictive power of the model is the key feature on 
which the reliability of the model results is predicated. 

e) turbulence models. These models are based on the Navier— Stokes equations and thus, at 
least in principle, have a built-in degree of generality, resiliency and predictive power. In 
particular, one-point closure turbulence models have the longest tradition (Chou, 1940, 
1945) and, naturally, they were the first tools used to construct the turbulence variables 
needed in O-GCM. The pioneering work in this field is that of Mellor and Yamada (1982, 
MY; Mellor, 1989) who applied the 1980 state-of-the-art turbulence modeling to 
geophysical problems. However, almost 20 years later, most of the work in ocean turbulence 
still employs essentially the same 1980 state of the art turbulence modeling (Rosati and 
Miyakoda 1988; Galperin et al.1989; Gaspar et al. 1990; Blanke and Delecluse 1992; Baum 
and Caponi, 1992; Ma et al. 1994; Kantha and Clayson 1994; Burchard and Baumert 1995; 
D'Alessio et al., 1998). This means that the basic difficulties that plagued the original MY 
model are still present. They can be identified with the closure of the pressure correlations 
terms, and the need to express the non-local third order moments (TOM) in terms of lower 
order moments. The pressure— correlations bring about adjustable constants which the 
model is unable to predict and the TOM are treated with the down gradient approximation 
which Moeng and Wyngaard (1989) have convincingly shown to be seriously inadequate 
(see paper I for details). 

f) Non-turbulence models. On the grounds that the above uncertainties have not been 
satisfactorily resolved, Large et al. (1994) have suggested an alternative model (called 
KPP), not based on any turbulence model. 
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g) New turbulence models. It is generally agreed that, potentially, turbulence models have 
far greater generality, resiliency and predictability power that any ad-hoc-model. In 
papers 1 and II, we offered the most updated 1— point closure models that one can construct 
today. Here, we offer a new 2 -point closure model that has been assessed over some ~80m 
turbulent statistics pertaining to large variety of flows. We derive the expressions for K m , 
and K g and run the 0— GCM with the new K's. The new model has no adjustable 
pa r ameters. 

II. New turbulence model 

The turbulence model to be employed here (Canuto et al., 1996—1999, cited as CD) 
was tested on about 100 turbulent statistics belonging to a wide variety of flows: pure shear 
flows, Rayleigh-Benard convection, rotating turbulence, 2D turbulence, rotating 
convection, freely decaying turbulence, stably stratified turbulence, etc. In all cases, the 
model reproduces the data quite accurately and thus the model can be considered 
validated. In addition, t he model has no adjustable parameters. Since this is the first 
application of this new turbulence model to ocean turbulence, it is necessary to begin by 
presenting the physical foundations of the model. 

Recent advances in phenomena related to turbulence (e.g., critical phenomena) have 
provided us with new tools and new perspectives on how to tackle the problem of 
turbulence. Like in many other fields, one begins with the case of homogeneous-isotropic 
turbulence for which an exact solution of the NSE found by Wyld (1961). The solution for 
the energy spectrum E(k) can be shown to be equivalent to the solution of a stochastic 
Langevin-type dynamic equation for the turbulent velocity fields u(k) of the form 

|,Uj(k) = ff l (k) + f‘(k) - ^ d (k)k 2 Uj(k) (4a) 

where E(k)~k 2 u(k)u(k). Eq.(4a) has an appealing physical interpretation: external forces 
are represented by the first term whose specific form is obtained from the NSE for arbitrary 
flows (shear, convection, rotation, etc). The second term is a random turbulent force 
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representing the action of the large eddies on a given eddy k even when the external forces 
are no longer operative (for example in the Kolmogorov region) and finally, the last term 
represents the dynamical viscosity felt by an eddy k due to all smaller eddies. Though 
derived for homogeneous and isotropic turbulence, Eq.(4a) seems to capture the main 
ingredients that one needs to build the canonical picture of a supply, transfer and enhanced 
dissipation of energy that is at the base of Kolmogorov original picture of turbulent 
dynamics. The challenge is to give an explicit expression to all the terms and prove that in 
spite of its "idealized" origin, Eq.(4a) is valid quite in general and capable of reproducing 
real flows. It must be stressed that we are not the first to suggest a Langevin-type equation 
as a starting point to describe turbulence. Several previous models were in fact recast in 
terms of a Langevin— type equations but, to the best of our knowledge, this is the first case 
in which Eq.(4a) has been shown to be a viable, successful, operational tool to treat real 
turbulent flows. To construct the explicit form of the dynamical viscosity, CD first 
generalized the well developed and well tested renormalization group method (RNG) to 
derive the following expression for ^(k) 

^(k) = + v = (v* + i / E(p)p- 2 dp)* (4b) 

k 

where u is the kinematic viscosity and ia (k) is the so-called turbulent viscosity. This 
expression tells two things: first, the integration from k to ® means that t'.(k) is 
contributed by all eddies smaller than k, as indeed expected on physical grounds and that 
while the Langevin equation is linear in the velocity fields, its "coefficients" like i/_.(k) are 
non linear since they depend on E(k). Thus, while in the NSE the coefficients are numbers 
and the velocity fields enter non-linearly (giving rise to closure problems), in the Langevin 
equation the situation is reversed. This a key point that allows us to derive a closed system 
of equations for the second— order moments like Reynolds stresses, etc. Next, we consider 
the first term in (4a). It represents the external forces driving turbulence, is a linear 
function of the turbulent fields and its form can be derived directly from the NSE. In the 
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present paper we consider the turbulent fields of velocity, temperature and salinity Uj, 9 
and cr in which case is the sum of three terms 

cO . c<7 (4c) 


ff f“ + {\ + f* 


where 


f“ = - ujjUjik) + 

,6 


f- = - a/ijtk) gjW 
<1 = “s P ij( k > fpl k > 


(4d) 

(4e) 

(4f) 


Here, U- • is the mean velocity gradient, and the a ' s are the expansion coefficients (T and 
bj 

S & are mean temperature and salinity) 


a T~~ pA % °S ~ pA % 

d 


i dp 


(4g) 


which must be computed using the equation of state for seawater. The superscript j. in (4d) 
means that only transverse components are allowed, that is, for any tensor A, 


A ir p > 


, P .(k) 
:m nap ’ 


(4h) 


where P.j=<f.j-kjkjk' 2 is the projection operator. Analogously, for the 9 and a fields, we 


have the Langevin equations: 

■^0(k,t) = f^(k,t) — y^(k,t)k 2 0(k,t) + f^(k,t) 

•^.<r(k,t) = f^.(k,t) — K^(k,t)k 2 a(k,t) + f^k,!) 

with 

X d (k) = X + X t 00, * d ( k ) = K + M k ) 

where x and « are the heat and salt kinematic diffusivities. The stirring forces are: 

f t ( /k,t) = -/jfu i (k), f‘(k,t) = - ^Uj(k) 

where the /?'s are the gradients of the mean temperature T and the mean salinity S a , 


„6_ cTT 
& 


(f _ 0s a 
Pi 


(5a) 

(5b) 

(6) 

(7a, b) 
(7c) 


1 £7 X . ’ ' 1 

As Eqs.(la-c) show, what one needs in practice the Reynolds stress r- which are the 
integral over all wavenumbers k of the two-point Reynolds stresses 

R j .(k)(5(k+k l ) = <u i (k,t)u j (k , ,t)> s (8a) 
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(8b) 


where s means symmetrization with respect to k and k' and i,j. The other required 


second-order moments are the temperature and salinity fluxes: 

jf(k)tf(k+k') = <u.(k,t)fl(k , ,t)> g (9a) 

j?(k)#k+k') = <uj(k,t)a(k’,t)> s (9b) 

As it will be shown, these variables entail other turbulence variables, and precisely, the 
variances of both the temperature and sait fields, 

e^k^k+k') = *<0(k)0(k')> ( 10a ) 

e (k)A(k+k') = i«T(k)a(k')>= (10b) 

as well as the temperature— salinity correlation 

e 0a (k)6{k+k') = <6{k)a(k')> (10c) 


The spectra of energy, temperature and salt fields E(k), ^(k) are obtained by 

integrating over the directions of the vector k (times a factor k 2 ) the density functions in 
the left hand sides of Eqs.(8)-(10), namely: 

E(k) = ,k 2 /R ij (k)dn k (11a) 

E <U (k) = k2/f W k)dn k (lib) 

The dynamic equations for the second— order moments (8)— (10) are derived by multiplying 
the stochastic dynamic equations (4a) and (5a, b) by Uj(k'), 0(k') and cr(k') and then 
averaging. The key difficulty is represented by the correlation of the first terms in (4a) and 
(5a, b) with the turbulent fields which physically represent the work performed by the 
turbulent forces f l (k). To compute these terms, we adopted two assumptions which are in 
agreement with the general understanding of turbulent processes: 1) energy generated at 
the largest scales is transferred locally to high k regions in a manner analogous to the flow 
of gas in a pipe, a model that applies to temperature and salinity as well, 2) as turbulent 
energy is transferred to smaller scales, the correlation functions of different fields, as well as 
correlation functions between different components of the velocity fields, decreases with 
increasing k even faster than the energy spectrum itself. The first assumption leads to a 
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representation of the energy flux n(k) of the form: 

n(k) = r(k)E(k) (12) 

where r(k) is the rapidity with which energy is propagated among different eddies. Eq.(12) 
can be viewed as the analog of j=pv with j-»II(k), p-»E(k) and v-*r(k). The energy flux fl(k) 


is related to the energy transfer function T(k) by: 

_ d TTM.\ 


T(k)=-g k n(k), /T(k)dk=0 (13a) 

T ( k ) = - r(k) ^E(k) - &r(k)]E(k) (13b) 

On the other hand, the function T(k) can be obtained by multiplying (4a) by u- (k* ) and 
separating the work of the non-linear interactions entering the first two terms in (4a), 

T(k) = A l (k)-2n t (k) k 2 E(k) (14) 

where the work of the force f* is defined by 

AJk) = k 2 /dn k dk'<u i (k',t)fj(k,t)> (15) 

Comparing (13b) and (14), we obtain: 

o k 

A t (k) = - r(k)| k E(k), *r(k) = /^(pjp^p (16) 

For the fields 6 and a, assumption 1) leads to relations analogous to Eqs.(15)— (16), namely: 

<tf(k')fj(k)> g = (47rk 2 )' 1 A^(k-fk') (17) 

«r(k')f^(k)> s = (47rk 2 )' 1 A ff ^(k+k l ) (18) 

where the expressions for the functions A^ (k) are analogous to Eq.(16), namely, 

A ^( k ) = - r «,.W|k E «,<T< k ) (19) 

? r ^(k) = / 0 X t (p)P 2 dp ( 20 ) 

*r a (k) = / 0 « t (p)P 2 dp ( 21 ) 

The functions satisfy the differential equation (F stands for x d and/or k^) 

f d “ T "d^d +F )-‘ (22) 

In the present paper, we shall be mostly interested in situations where 

v$>v, x d >X, * d > K ( 23 ) 

and thus it follows from Eq.(22) that 
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( 24 ) 


(*d- K d> = a t v d M 

where the turbulent Prandtl number <7 = 0.72 and ^(k) is given by Eq.(4b). As for 
assumption 2) above, it leads to the following relations: 

<f|(k,t)Uj(k',t)> = - (8irk 2 )' 1 A t (k)P j j(k)£(k+k') 

<u i (k')«(k')f* ftff (k)> = 0 
<«(k’)f! iff (k)> = 0 

<cr(k')f| s ik)> = 0 (25a) 

Finally, the dynamic equations for E(k) is: 

ItE(k) = A ext (k) + A,(k) - 2k 2 i' d (k)E(k) (25b) 

The information we have derived concerning the terms in the basic equations (4a) and 
(5a, b) is sufficient to allow us to derive a closed system of equations for the second-order 
moments (8) — (10). 

III. Dynamic Equations for the Second— Order Moments 

Multiplying the stochastic equations (4a) and (5a, b) by the fields Uj, 0 and a, it is 
straightforward to derive the following dynamic equations for the second-order correlation 
functions: 

*!t R ij(k) = ( 8 ”k 2 )- 1 A t (k)P ij (k)-[St m R mj (k)] s - 
- ![V X , RfkJljj + * [DRtkJIlj - kV d (k)R..(k) 

- a T (gJj^( k )] s + a s (g|jj( k )] s (26) 

+ p/dOlf 

-k> d (k) + x d (k)]jf(k)-^R im (k) 

[-2o T e^ k ) + tt g e <r^ k )] g j P ij^ k ) ( 27 ) 

|^k) = (4ffk 2 )‘ l A d (k) - j3fjf(k) - 2k 2 x d (k)e S (k) (28) 

1 e^(k) = - y9fjf(k) - 0?jf( k) - 2k 2 [x d (k)+ Kd (k)]e^(k) (29) 
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In addition there are two equations which can be obtained from Eqs.(27) and (28) with the 
substitutions: 


jf^jp V _a s’ W V A <r’ Vd 

In Eqs.(26)-(27) we have introduced the tensors 

u. .=S..+V.., 2S..=U. .+U. 2V..=U. .-u.. 

i,J ij ij ij i,J J,i’ J] bJ Jh 

where S- j and V- are the mean shear and vorticity and D is the operator 

D = k j( Sp Vj^ 


(30) 

(31) 

(32) 


Eqs.(26)— (29) must be solved together with the incompressibility condition 

k.R..(k)=k.R..(k)=k./ a =0 (33) 

i ij v ’ j ij v ; i J i ' ' 

This suggests that one may solve the above equations in spherical coordinates where one 

can write (for any j) 


j(k) = ykfy + j^kfy (34) 

where e^ ^ are the basic units vectors. Due to condition (33), the only non-zero 
components of Rjj lie in the plane perpendicular to k so that their number is three and the 
general decomposition has the form: 

R y (k) = e(k)P jj (k) + a(k) Pjj (k) + b(k), sj (k) (35) 

where e(k) is the energy density in k-space while the orthogonal tensors p and rj should be 
chosen so as to avoid singularities in the functions a(k) and b(k) which invariably arise 
when the tensors p and tj depend on 0 as #-*0. This point and the details of the 
decomposition are fully discussed in Sec. Ill of Canuto and Dubovikov (1996b). In the 3D 
case, the final equations for the scalar functions e(k), a(k), b(k) etc are quite cumbersome. 
Thus, we shall consider only the case of direct physical interest: the mean fields U-, T, and 
salt depend only on the coordinate where there is a non-^zero gravity. The final result is: 
^e(k) = (4^k 2 )‘ 1 A l (k) - 2k 2 i' ( j(k)e(k) + M^k) + M a a(k) + M^b(k) - 

- 4£ (k) “ A W k) + A ^ (k) + A K (k) (36) 
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(37a) 


f^k) = - 2 k 2 i' d (k)a(k) + M fl a(k) + M a e(k) + M ab b(k) - 

+ + V$ (k) + 6 $ k) 

with 

<5 = A^cos 20 +A^sin 2 <^, - A^cos 20 +A^sin 2 <^ 

S^= - A^cos 2 ^>-A^sin 2 <f>, b = A^cos 2 <^>— A^sin 20 (37b) 

|b(l) = - 2k 2 i^ d (k)b(k) + M o b(k) + M b c(k) - M ab a(k) - 

+ V*(k) + l,j^) + 7 3 jJ(k) + 7,j>) (38a) 

with 

7 ^= A^sin20-A^cos20, 7 ^= - A^sin 2 </>-A^cos 20 
7 ^= - A^sin 20 +A^cos 20 , 7 ^= A^sin 2 <^+A^cos 2 </> (38b) 

f|jjj(k) = - 2kV d (k)i^(k) +M^-Sz(l-sin2^)j^-/3^(k) - 

- 2 A^(k) + A^k) + ( a(k) + ( ,b(k) (39a) 

^ = /?^sin20-/?jjcos20 , = - /?^sin20-/?^cos20 (39b) 

m4 ( k) = ‘ 2k2 P t 'd< k )4 (k) + M <4 + Szi o~ ^ e(k) “ 

- 2 A^k) + A Je^k) + , a(k) + , 2 b(k) (40a) 

^ = /^cos20+/?^3in2 0 ,7 ^ - /i^cos20-t-/?^sin20 (40b) 

5 ,e 9 !k) = (4irk 2 )- 1 A^k)-2k 2 <rj 1 i/ d (k)e^(k) 

(«) 

a'W 11 * = - " ^ k > ~ 

- 4 ?( k > -/$*«-(#<) <«) 

We have defined the following variables (z=cos$): 

= iS(l— z 2 )cos20 + D (43) 
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(44a) 


-4S-'1V1 = (1+z) 2 + (1 -z) 2 cos4$ 

a 

— 4S _1 M^ = ( 1— z) 2 sin4 0 
S J M b = (1— z)(l-sin20) -z 
= — £Sz 2 cos20 + D, M^= -jScos20 + D 
2S _1 D = cos2 4> ( 1 — z 2 ) (k^ k -z^ z ) + (l-sin2^)^ 

p = 4(l + a^) =1.2 (44b) 

/?^ ,<7 = - 2'^(sin<?!>-cos0)/?^ ,<T 

(3^ a = 2 _ ’z(cos<^ + sm<j))0^ ,(J (44c) 

\^° = - 2'*z(sin</»+cos^)ga T g 

A*' = 2'*(sin$-cos0)ga T g (44d) 

Finally, the equations for 

jj(k), jj(k), e,(k) (44e) 

are obtained from Eqs.(40)— (41) via the substitution (30). In the above and subsequent 
equations, 

s = PSyS,/’ N2 = go i^’ Ri = ^ (44f) 

S, N and Ri are the mean shear,the Brunt-Vaisala frequency and the Richardson number. 


IV. Spectra of the Second— Order Moments 

The above equations do not contain free parameters and can be solved numerically 
once the initial density functions are known. Concrete examples were discussed in Canute 
and Dubovikov (1996c) for different types of shear driven flows for which we can compare 
the model results against DNS/LES and laboratory data. The model reproduces such data 
quite well. The case of buoyancy driven flows was also studied (Canuto and Dubovikov, 
1997a) and in that case too, the model results were in good agreement with a variety of 
data. This performance of the model gives us confidence that the extension discussed here is 
based on solid foundations. Our final goal is to express the one— point second order 
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moments like the Reynolds stresses, the temperature and salinity fluxes in terms of the 
gradients of the large scale fields. Since a numerical solution of the above dynamic 
equations would not be of practical use, especially if the results are to be used in an 
O-GCM, an analytical model is needed which by necessity requires that we make some 
approximations. We shall employ the same philosophy that was employed earlier when we 
treated shear and buoyancy driven flows. The gist of the method is as follows. In the 
inertial ranges, Eqs.(36)— (42) can be solved perturbatively since we have the small 
parameters 

S(k 2 i/ t ) _1 , gft T /^(k 2 * t )- 2 , ga s p a (k 2 K t )- 2 (45) 

Once the inertial spectra are obtained, we extrapolate them all the way back to the 
infrared cut-off k^ below which we assume that E(k)=0. Upon integrating the resulting 
spectra over all k's, we obtain the desired one-point correlation functions. In the case of 
shear— driven flows, the results obtained with this method differ less than 15% from those 
obtained by numerically solving the full equations (Canuto et al., 1999c). We shall carry 
out the procedure without salinity which we shall incorporate in the next paper. To 
calculate the spectra, we first need to integrate the k— space densities, Eqs.(8)-(10), over all 
directions of k. In turn, the densities are expressed in terms of the variables in Eq.(36)— (42) 
as described in detail in Canuto and Dubovikov (1998b). Furthermore, we shall consider 
first the stationary case. To the zeroth-order in the smallness parameters (45), the only 
non— zero results are the isotropic components of the energy densities e(k), e^(k) and e^(k) 
which are obtained by retaining only the first two terms in Eqs.(36) and (41) and in the 
equation for e^ which is obtained from (41) and the substitution (30). The results are: 


47rk 2 e(k) = E(k) = Kof 2 / 3 k- 5 / 3 

(46) 

4^k 2 e^(k) = E^(k) = Bac^r ^ 3 k' ^ 3 

(47a) 

47rk 2 e (k) = E (k) = Bacr^k' 5 / 3 

(47b) 

Ko = lj, Ba 

(48) 
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where Ba is the Batchelor constant and the c's represent the rate of dissipation of the total 
kinetic and potential energies for the temperature and salinity fields. To compute the first 
order terms in the smallness parameters (45), we must consider the first and third terms in 
(37a), (38a), the first, fourth and fifth terms in (39a), (40a) and in the corresponding 
equation for j , whereas in (36) the second and third terms should be taken into account 
provided that only the component ene-e is retained in the second term since the e-part has 
already been used in the zeroth-order. The results of the above procedure, as well as 
further intermediate steps, are presented for the simplified case in the absence of salinity, 
whereas the final results will include salinity. 

Reynolds stress spectrum 

R(k) = ^IE(k) + C a (k)a + + C (k)t + C A (k)A (49a) 

w r here the basis tensors are: 

a = S^S, £ = S- 2 S 2 - 

t = S- 2 [V,S], A = g' 2 gjgj- jl (49b) 


The coefficients C's are given by (superscripts indicate the perturbative order at which 

0 

they were obtained and 0=/3 )): 


C;(k)=-^Sr ! /3 k - 7 /3 

(50a) 

C 3 (k) = 0.92(1 - a^r^g/T 1 )^- 1 / 3 ^ 11 / 3 

(50b) 

C t 2 (k) = 27 S 2 k- 3 

(50c) 

Cj*(k) = -1.08 (1 - a t c^€' 1 a T g/?' 1 )N 2 S 2 c' 2 / 3 k' 13 / 3 

(50d) 

C|,(k) = -3.3S 2 k- 3 

(50e) 

C 2 x {k) = -1.3N 2 (1 - 1.44e^r 1 ft T g/?- 1 )k- 3 

(50f) 

where the Brunt-Vaisala frequency N is: 

(Mg) 

N 2 - 

(50h) 

Heat flux spectrum: 


J(k) = J (k)g + JU 

w g v /b 0 U 0 

(51a) 

g=8-'g, U = U-‘U 

(51b) 
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(52c) 

(52a) 


J‘(k) = 0.93(3 (1 - 1.44^r 1 a T g^ 1 )t^ 3 k' 7 ^ 3 
J u (k) = ^P' 1 /^ + S P 4 0 " 1.44^r 1 Q T g^- 1 )]Sk- 3 

V. One-point Reynolds stresses and temperature flux 

The case of the Reynolds stresses in the presence of shear was discussed in detail in 
Canuto and Dubovikov (1999c) and we shall not repeat it here. The final result for the 
one-point Reynolds stress, arrived at by integrating (49a) over wavenumbers from the 
infrared cutoff to infinity, is: 

T if ¥j-i K4 ij ,53a) 

r . = C a + C.t + C x A (53b) 

lj G 1 A 

where the basis tensors were defined in Eqs.(49b). The coefficients C's are : 

C <r = -5 StV V /3<r m (54a) 

Ct-STSV'm (54b> 

C A = - 0.65N 2 k Q 2 (l - 1.44e^- 1 a x g/?- 1 ) (54c) 

a m = 1 - 0.62(1 - 0.72f^r 1 o T g/?- 1 )N 2 r 2 / 3 k o 4 / 3 (54d) 

Upon integration of the thermal spectra (51) and (52), we obtain the two components: 

J g = 0.7/?f 1 / 3 k o 4 / 3 (l - IMetf-'aygP- 1 ) (55a) 

J u = 0.156 /3Sk* 2 [l + 1.43(1 - (55b) 

We recall that k is obtained via the relations: 
o 

K = /( E(k) dk = | (55c) 

k = (Zf/hK- 3 ! 1 (55d) 

0 * 

Some modifications are however required. As stated earlier, the above results correspond to 
the stationary case and thus, they imply that production=dissipation, that is, 

P=e, p 0 =f 0 ( 56 ) 

where P and represent the rate production of turbulent kinetic energy and of 
temperature variance (that is, potential energy). To introduce the P/e dependence into the 


22 



equations, we must somehow keep the left hand side of Eq.(36). In carrying out the 
perturbative approach, we obtain the following differential equation for the functions e,a,b 
to the first order in the smallness parameter (45): 

| t f(k)= £ ‘/3 k 2 /3[ f « (k) _f (k) j (57) 

where f^ is the first order term representing any of the functions e,a,b in the stationary 
limit. The same equation holds true for the first order of the spectrum C^k). To integrate 
(57) over k for the component C^.(k), we assume that its form is the same as in the 

7 / 

stationary limit, that is, k" /3 . Then, from (57) we obtain: 

♦'• 1/3k „ 2/3 lt c i< k ) = c J- c i( k ) < 58a > 

Following a practice widely used in treating the Reynolds stresses, we shall assume that 


ft c ;< k > = c o kk -' - cpvtp-o 


Thus, Eq.(58a) gives finally 


(58b) 


Cgk) = FCj(k) (58c) 

F-> = 1 + i(P/(-l) (58d) 

Proceeding in an analogous fashion, we obtain 

J g = J g p , (59a) 

P = 1 + 0.084 [P/c-1 + ^(P^/e^-l)} (59b) 

We can continue the analysis of the P/f dependence to the next order of perturbation, as it 

was described in detail in Canuto and Dubovikov (1999c). The results will be presented 

below but before we must discuss another correction to C 1 which is due to the fact that an 

o 

unmodified extension of the initial expression for the first-order function e^k) up to k^ 
results in a violation of the positiveness of the energy density e(k) in the region k o <k<k^ 
where 


k = 


i 



x 


KS 


Such a violation can be prevented by multiplying the function e^(k) by (k/k^) 
region k o <k<k^. This leads to the following correction factor to C^(k) 


2 / 3 


(60a) 
in the 
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(60b) 


F c = ^ - 3x _1 (l - ygX' 1 ) for x> 15/8 
F c =l, for x<15/8 

Next, we substitute the variables from Eq.(55d) and adopt P q to obtain 

(#= 0J g (60c) 

Substituting (55a), we obtain a linear equation in e^/e whose solution is 

-f = y^ 0 /?N- 2 x 2 Ri(l + O.lGx^i)- 1 (61) 

The last modification we have to carry out accounts for corrections of higher orders in 
(N 2 ) n to the function o > Eq.(54d) where we took into account only the first correction 
and thus a has the form 

%= W (62a) 

Since a m >0. a consistent way to extrapolate (62a) to the regions where £~1 in the spirit of 
Pade 1 approximations, we write 

l-e-O+O' 1 (62b) 

VI. Reynolds stresses, temperature and salinity fluxes 


After performing the modifications indicated above, we obtain the following results: 

Reynolds Stresses: 


T if u i u rH 

(63a) 

r ij = C o° + C t l + C A A 

(63b) 

where the basis tensors are defined in Eqs.(49b) and: 


C „ = - 4 f SS m 

(63c) 

C. = 47 ^ rS 2 S 
t '3 c m 

(63d) 

C, = - 8 3-^7 — 2 fN 2 S, 
A 4 c h 

(63e) 

t=2Kc 1 

(63f) 

In particular: 


— v dU — v dV 

uv, = - K mW' vw = " I W 

(63g) 
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Vertical Heat Flux J . 


K = 2 - 2 S 

m f m 


J 6 -Wd- - ~K h ^ 
K h = 2 7 S h 


Vertical Salt Flux J : 

<7 

K =2- 2 S 
a to 

Dimensionless Structure Function for Momentum: 


where 


S m = 25 F , Vm 


P c= 5 F c = 1 ^ x l+ t6 x ' 2 for 
P c = g, for x< 15/8 


where the dimensionless variable x is given by 


The other dimensionless functions are: 


a m = 1 + ( W F c ! 

1 0 2 77^ = 5f' 2 F 3 X 2 Ri[l+ Y2^ 1 F’ 1 (5+7F 2 P^ , )(l-0.16x 2 f : ' j R^Ri)^] 
10 2 r j(j = - 5f' 2 F 3 X 2 R /) Ri[l+ i 2 P i Fj 1 (5+7F 2 f' 2 1 )(l+O.16x 2 P i Ri)0 i ] 


0- 1 = 1 + 0.16 f^x 2 (l— R )Ri - 0.013 F 2 x 4 R^Ri 2 ^ 
<j) A = 1+ 0.08 x 2 (l— R^)Ri 

^=(6 a sl 


Furthermore: 


S I rsj 4 rv rv 

m = 35%' % = 1+ V % 

10 2 ^ = 4P 3 F 4 X 2 Ei[l+ { 2 F 2 F-'(5+7 F j P 3 i )(1-0.16x 2 P R^Ri)^] 
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10 2 ^ = - 4F 3 F 4 x 2 Riiyi+ ^ 2 $F-J(5+7FPp(l+0.l6x 2 $m)il>J (67c) 

t 1 = 1 + 0.16P x 2 (l— R )Ri - 0.013 P 2 x 4 R Ri 2 tf> (67d) 

2 2 P 2 P 

The functions 7's and F n are given by: 

7 = 0.12F F , 7 = 0.095F (68a) 

3 1 2 4 2 

F n > = l + fn(l+n)-‘(f-l) (68b) 

P n =1 + S n ( 1+n )' 1 (7- 1 ) (68c) 

Dimensionless Structure Function for T emperature and Salinity: 

S h = O.OSeF^Pj 1 - 0.08x 2 ^RiRp^ (68d) 

= 0.056P 2 (P 4 + 0.08x 2 </>Ri)^ (68e) 


VII. Local and Non-Local Models 

The above expressions require the knowledge of both the turbulent kinetic energy K 
and of its rate of dissipation c. In principle, K and c satisfy the non-local dynamic 
equations: 

w + i F(K > = - T ij S ij + f < 69a > 

Ir + lz F (') = l _c 3 r ij s ij + c /* a T J » _ga s J <7M tK " ,- V 2K 1 

(69b) 

The values of the coefficients c have been discussed at length by Burchard and Baumert 

1 > 2^3 

(1995). F(K) and F(c) denote the fluxes of K and t and they are third-order moments. The 
most complete expressions for these third-order moments (TOM) are those given by 
Canuto et al. (1994) which were successfully tested in the PBL (planetary boundary layer) 
against results from LES (large eddy simulation). Use of these TOM, together with the 
above model for the diffusivities, would represent a very complete treatment of ocean 
mixed layer turbulence. It is the K-e model widely used in engineering turbulent flows. 

Here, we present numerical results for the local limits of (69a, b) which we can write as 

- + gmJ £ = c (69c) 

c = K 3 / 2 A 4 (69d) 
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where A is a mixing length. The first equation has the simple physical interpretation of 
production=dissipation, while (69d) is a statement of Kolmogorov law. Using the above 
expressions for r- and J^, Eq.(69c) becomes the equation for the dimensionless variable x 
defined in Eq.(66c). We obtain: 

X* - *(s m - His,,)-* (69e) 

where both S m ^ depend on x. The Richardson number Ri is defined in Eq.(44f). For a 
given Ri, x is thus uniquely determined, x=x(Ri) and so are the structure functions 


. = S_ , (Ri) 
m,h m,h v ' 

If we use the other local relation (69d), we finally obtain 


(690 


K . = 2A 2 Sx-’S . 
m,h m,h 


(69g) 


which yields the turbulent diffusivities in units of A 2 S. Finally, we note in the 
production=dissipation local case, there is further simplification: 


F =P =1 (69h) 

1 *? 2 ? 3 ? 4 1 ■? 2 ^ 3 

We must stress, however, that the use of a local model is not required by the turbulence 

model we have presented being made here only for reasons of simplicity. The length scale A 

is determined using the Deardorff-Blackadar formula: 

A = 2- 3 / 2 B C, B =24.7, 

1 1 

l — min(ij|, £ ) (69i) 

l = kz l ((. +/tz) 4 , l = 0.17H (69j) 

10 0 0 

where 4q 2 =K is the turbulent kinetic energy, N is the Brunt— Vaisala frequency, «=0.4 is 
the von Karman constant and H is the mixed layer depth. When used within the NCAR 
CSM Ocean Model, H is determined as the depth where the buoyancy difference 

g[p(H)— /9(surface)]p(H)' 1 = 3 10' 4 ms' 2 (69k) 


VIII. Ocean GCM 

We tested the new vertical diffusivities in a global ocean general circulation model, 
O-GCM. We used the NCAR CSM Ocean Model produced by the University Corporation 
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for Research. National Center for Atmospheric Research, Climate and Global Dynamics 
Division. They developed their model by modifying the MOM 1.1 GFDL code as described 
in the NCAR CSM Ocean Model Technical Note, "The NCAR CSM Ocean Model", by the 
NCAR Oceanography section. We employed the stand-alone 3°x3° configuration of the 
model detailed in their technical note with the default parameter values. It has 3.6° spacing 
in longitude and a variable spacing in latitude increasing from 1.8° at the equator to 3.4° at 
17° N. S and then decreasing back to 1.8° for 60° N, S and poleward. There are 25 levels of 
increasing thickness in the vertical, with the surface level 6 meters thick. The option for the 
GM mesoscale eddy parameterization was enabled. Bulk forcing with a seasonal cycle plus 
a 1/2 year timescale restoring condition on the salinity is used, except under sea-ice where 
there is strong restoring. This configuration corresponds to the case B— K described in Large 
et. al. (1997). It should be noted, however, that for determination of the length scale in our 
turbulence model we used the program's definition for mixed layer depth (a buoyancy 
difference from the surface of 3 10‘ 4 ms' 2 ), which is different from that graphed as a 
diagnostic in Fig.5 of Large et al.(1997). We initialized our runs with annually averaged 
Levitus data and ran for 126 momentum years. As in Large et al. (1997) a 3504sec timestep 
for momentum is used, while for the first 96 momentum years the tracers are accelerated 
by a factor increasing from 10 at the surface to 100 for the deep ocean. We then set all 
timesteps equal for the remaining 30 years as they did. 

First, we ran the NCAR program as is, with the option for the KPP mixing enabled, 
producing the KPP data presented in the figures below. Then, in place of the KPP module, 
we inserted a module which uses our new model for the diffusivities for momentum and 
heat with the salt diffusivity set equal to that of heat. To save computing time, we 
constructed tables of the dimensionless functions S m ^ and of the dimensionless variable y 
(obtained from solving Eq.69e), 



vs. Ri. Then, for each point in space and time these were interpolated to the local Ri. To 
construct the diffusivities K m ^(model), we used the results (63h) and (64b) which can be 
rewritten in terms of (70a) as 


K 


JfS = *B y*S m . 
m,h y 5 \ J m,h 


(70b) 


IX. Below the Mixed Layer 

Below the ocean mixed layer, the external wind— generated shear is too small to 

generate turbulent mixing and yet, even in regions where both the temperature and the 

salinity gradients are stably stratified, it is usual to assume "background diffusivities" for 

viscosity, heat and salt diffusivity which are believed to be caused by internal wave 

breaking (Large et al., 1997). It would be preferable not to do so but rather model the 

physical processes causing this background mixing. Our main assumption is that the 

turbulence model has given us the correct functional dependence of the K m b g on Ri and 

R and that such diffusivities can thus be used below the ML. Since all the arguments 
P 

discussed below, are valid for any of the three K's, we shall use only the generic symbol K 
and write succinctly 

K = K(Ri,R ) (71a) 

The key problem is how to define and thus compute Ri. Here, we shall make us of the 
measured data (Gargett et al., 1981) on the vertical shear generated by the wave breaking 
phenomenon. By integrating over all wavenumbers one can compute the shear due to 
internal waves, S,. One can then form a corresponding Ri wb as follows: 

Ri wb = N 2 /S^ b (71b) 

where 

N 2 = -g^ | (71c) 

Gargett et. al. (1981, sec. 5) confirmed earlier arguments by Munk (1966) that Ri wb ~l. To 
those argument, we would like to add the following consideration. As the value of Ri cr , 
above which there is no longer turbulent mixing, computed from our model is 0(1), if Ri wb 
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were >>1, there would be no turbulence generated by the internal waves at all. On the 
other hand if Ri w ^ were <<1, there would be a very strong turbulence producing a 
viscosity sufficient to damp out the waves themselves. The wave— generated turbulence is 
thus self-limiting. Since the turbulence model gives a precise value for Ri cr , while the 
above argument only tells us that Ri w | 3 ~0(l), we shall write: 

Ri , = c Ri (71d) 

wb cr 

where c is a constant reasonably close to unity. We have found that c=0.88 gives a 
diffusivity close that measured by Ledwell et al. (1993). Since in the local model, the K's 
are also proportional to the length scale A or see Eq.(69i,j). Below the mixed layer, we 
thus need an analogous 1 ^. We shall use the same formal expressions (69i,j) but with 
different / ( wb) which we compute as follows. Assuming a Kolmogorov spectrum at 
wavenumbers upward of a breakpoint and integrating, we obtain: 

* (wb) = (3Ko) 3 / 2 (B k o ) 4 (71e) 

where Ko=1.6 is the Kolmogorov constant. We identify k^ with the best value of Gargett 
et al. (1981) for the break in slope of the observed spectrum of internal waves, namely 

k^ = yq 2 7r radians/meter (71f) 

Thus, f Q (wb) is knowm and so is Similarly, y w ^ is obtained by solving the 

production =dissipation, Eq.(69c). Thus, the complete wave-breaking expressions for the 
three diffusivities are: 

K m,h,s< wb ) = * B /( wb ) S wb y wb S m,h,s< Ri wb’ V (72a) 

We add together the diffusivities calculated using the shear resolved in the ocean model 

and the background diffusivities, ensuring continuity in the transition between regions 

where external excited shear dominates and those where the internal wave shear does. We 

thus take the total diffusivities to be: 

Vh,s= W a ’V + W wb -V (72b) 

In the statically unstable case (Ri< 0), we set K m ^ g (wb)=0. The very large mixing due 
to convective instability makes the background irrelevant in this situation in any case. 
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X. O-GCM results 

In Fig.l we present the dimensionless structure functions S m ^ vs.Ri. They are given 
by Eqs.(66a) and (68d). The graphs correspond to the case without salinity. In Fig.2 we 
present the ratio K m /K^ (again for the case without salinity) versus the LES data (dots) of 
Wang et al. (1996). In Figs. 3-8 we present various diffusivities and their ratios for the full 
case with salinity and shear. The latter is represented by Ri, defined as the ratio of N 2 , 
Eq.(50h) and the shear squared. The numbers on the curves represent the Turner number 
R defined in Eq.(66i). The nomenclature SF, DC, DS and DU was introduced and 
discussed in the Introduction. 

In Figs. 9-20 we present the results of the O-GCM which comprise both global 
temperature and salinity profiles as well as for the case of different ocean basins. In 
Figs. 21— 29 we present the diffusivities in different ocean basins. Particularly interesting is 
the case of the Canary Islands since one can compare the model results with the value of 
0.11±0.02cm 2 s- 1 measured by Ledwell et al. (1993). In each case, we compare our results 
with Levitus data as well as with the results we have obtained by running the same code 
with the KPP model which is available only for the case K g =Kj i - 

Finally, in Fig. 30 we present the polar heat transport. 

XI. Discussion 

It seems natural to require that the reliability of a model used to describe a given 
phenomenon should precede rather than follow the application of the model. Stated 
somewhat differently, a model should be judged on more than its performance in a given 
instance. Regrettably, the 1—point closure models used thus far to treat vertical mixing and 
following the Mellor— Yamada seminal work, have suffered and still do, from the presence of 
several adjustable parameters which have not allowed a clear appraisal of their performance 
before their use in the ocean context. In that respect, the new 2-point model presented here 
is quite different for it is being used only after it has exhibited a pedigree based on its 
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performance on a wide variety of turbulent flows. Its credibility, manageability and 
resiliency makes the best candidate for the computation of the vertical diffusi vities. 

The work is however not complete. In fact, it is known that double— diffusion 
phenomena change a smooth gradient into one that exhibit a step-like structure. As of 
today no theory that we know of is capable of encompassing such phenomenon. It is a 
challenge that we shall undertake in future work. 
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Figure caption 

Fig.l The momentum diffusivity K m vs. Ri for different values of the Turner number 
The label DC and SF are defined in the Introduction. 

Fig. 2 Same as in Fig.l for the DU and DS cases. 

Fig.3. Heat diffusivity vs. Ri for different R^ for the DC and SF cases. 

Fig. 4 Same as in Fig.3 for the DU and DS cases. 

Fig. 5 Salt diffusivity K vs. Ri for the SF and DC cases 

o 

Fig. 6 Same as in Fig. 5 for the DU and DS cases 

Fig. 7 The turbulent Prandtl number vs. Ri for different R . The heavy line 

corresponds to the case of laboratory (see paper I, Figs. 3-^1). Cases DC and SF. 

Fig. 8. The ratio of K /K, vs. Ri for different R . Cases DU and DS. 

° m' h p 

Fig. 9 The ratio of K m /K g vs. Ri for different R^. DC and SF cases. 

Fig. 10 The ratio of K m /K g vs. Ri for different R^. DU and DS cases. 

Fig. 11 The ratio of K^/K g vs. Ri for different R . DC and SF cases. 

Fig. 12 The ratio of K^/Kg vs. Ri for different R . DU and DS cases. 

Fig. 13 The efficiency ratio T (paper II, Eq.59f) vs. Ri for different R^. DC and SF cases 

Fig. 14 The global ocean temperature using the O-GCM discussed in VIII-IX with 
the new diffusivities. The Levitus (1994) data are the solid line. We have also run the 
O-GCM code with the KPP model (Kg=K^) and the results are indicated by 
diamonds. The results with our new 1— point closure model with K g =K^ (paper II) are 
shown by squares while the present model results with are indicated by 

asterisks. 

Fig. 15. Global salinity 

Fig. 16. Artie ocean temperature 

Fig. 17. Artie ocean salinity 

Fig. 18. Atlantic ocean temperature 

Fig. 19. Artie ocean salinity 
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Fig. 20. Pacific ocean temperature 
Fig. 21. Pacific ocean salinity 
Fig. 22. Indian ocean temperature 
Fig. 23. Indian ocean salinity 
Fig. 24. Southern ocean temperature 
Fig. 25. Southern ocean salinity 

Fig. 26. The four diffusivities K , (cmV 1 ) for the Papa station 
Fig. 27 The four diffusivities K , (cmV 1 ) for the Artie ocean 

Fig. 28 The four diffusivities K , (cmV 1 ) for the Canary Islands. 

Fig. 29 Same as in Fig. 26 but for the first 1km 
Fig. 30. Same as in Fig. 27 for the first 1km 

Fig. 31 Same as in Fig. 28 for the first 1km. Ledwell et al. (1993) value of 0.11±0.02 cmV 1 
(see, however, the discussion in the main text) 

Fig. 32. Same as in Fig. 26 for the first 40m 
Fig. 28 Same as in Fig. 22 for the first 40m 
Fig.29 Same as in Fig. 23 for the first 60m. 

Fig. 30 Polar heat transport vs. latitude for three different models. 
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